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As computational fluid dynamics techniques and tools become widely accepted for real- 
world practice today, it is intriguing to ask: what areas can it be utilized to its potential in 
the future. Some promising areas include design optimization and exploration of fluid 
dynamics phenomena (the concept of numerical wind tunnel), in which both have the 
common feature where some parameters are varied repeatedly and the computation can be 
costly. We are especially interested in the need for an accurate and efficient approach for 
handling these applications: (1) capturing complex nonlinear dynamics inherent in a system 
under consideration and (2) versatility (robustness) to encompass a range of parametric 
variations. In our previous paper, we proposed to use first-order Taylor expansion collected 
at numerous sampling points along a trajectory and assembled together via nonlinear 
weighting functions. The validity and performance of this approach was demonstrated for a 
number of problems with a vastly different input functions. In this study, we are especially 
interested in enhancing the method’s accuracy; we extend it to include the second-orer 
Taylor expansion, which however requires a complicated evaluation of Hessian matrices for 
a system of equations, like in fluid dynamics. We propose a method to avoid these Hessian 
matrices, while maintaining the accuracy. Results based on the method are presented to 
confirm its validity. 


I. Introduction 

T O make a CFD solution practically and broadly useful, it must be fast, easy to set up, robust and accurate. 

Examples such as design optimization and exploration of flow scenarios with different inputs/parameters 
require repetitive calculations of solutions varying over the parameters space. Hence, it is enormously important to 
have a predictive model that possesses the following property: (1) accuracy, (2) efficiency, and (3) robust 
applicability. The first property demands that the model retains the accuracy of the full solution (whether it be Euler 
or Navier-Stokes equations based); the second property is to require only a fraction of the cost of the full solution; 
and the third property is to require that the model be applicable over a sufficiently large variation of parameters, 
even revealing some physical phenomena not seen in the baseline solution. 

The goal of the study is to take the traditional CFD goal to a level higher — beyond computing a solution 
corresponding to each set of parameters each time, but producing a model that can be used over different sets of 
parameters, varying over a sufficiently large space, many times. Simply put, the modeling is based on a given 
baseline trajectory (configuration), then the model so developed should be able to perform outside of this baseline, 
in a new configuration substantially different from the baseline. 

A key element for the modeling to be successful is to be able to preserve nonlinearity that exists in problems of 
interest, thus, requiring at least a nonlinear formulation in the modeling from the outset. The modeling principle 
should be applicable to a range of governing equations; for fluid dynamics they can be Euler or Navier-Stokes 
(including RANS, EES, and DNS) equations. 

The objective of this paper is to present the concept of our nonlinear modeling approach; we shall describe the 
basic concept and steps for constructing such a model and then demonstrate its efficacy of the above three properties 
for problems of practical interest in aerodynamics. Cases included in this paper are based on the Euler equations 
since our interest primarily lies in unsteady flow phenomena in which the main physical driving force is the inviscid 
one. However, the same approach can be applied to a steady flow problem where the input parameter can be a new 
stationary Mach number or a new flight angle of attack. 
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Let us consider a time-dependent nonlinear system, 


3x 

¥ 


= f(x,u) 


where X is a «xl state vector and U a mxl input vector. A common approach for constructing a nonlinear model 
to meet the above three properties for accuracy, efficiency, and robustness (versatility) is to formulate the problem in 
a reduced space in order to reduce the computation effort. An appropriate projection can be chosen to construct a 
nonlinear ROM in a reduced basis <j), expressed as 

X, =0^f(0X^,u) (3) 


where X^ is a rx\ state vector, with r « n. In Eq.(3), it is critically important to have an accurate nonlinear full- 
order model X before applying projection; a proper construct of the projection matrix ^ is a key step that affects the 
success of the modeling. Rewienski and White^*^ introduced the trajectory piecewise-linear (TPWL) approach for 
constructing a nonlinear model: including the projections of state variables and the governing nonlinear differential 
equations onto a subspace together with a proper weighting of a number of truncated linear models that are 
expanded about some operating points (or we preferably call sampling points in this paper). Gu and Roychowdhury^^^ 
proposed a similar trajectory-based approach by projecting the full order model into a low order nonlinear manifold 
subspace and demonstrated improved accuracy by their method for analog circuits and a biochemical system. 
Gratton and Willcox'^^^ combined the TPWL and proper orthogonal decomposition (POD) methods to derive a 
reduced model to study a flow control scheme for a supersonic diffuser. He et al.'^"^^ pointed out instability arising 
from applying the TPWL method with the POD method, or to over-fitting problems. They discussed that instability 
appear when including basis vector that corresponds to smaller eigenvalue in the POD method. 

Preserving nonlinearity of the problem under study is utmost important for obtaining a model’s predictive 
capability. The work by Rewienski and White ^*'on micromachine devices inspired us to use piecewise linear local 
solutions for study of nonlinear unsteady aerodynamics. In the context of control theory, gain scheduling is a very 
common practice to approximate a nonlinear control plant^^\ and it shares the same concept with TPWL, that is, 
divide the complex nonlinear model into sub-models^®l Computational cost can be reduced by gain scheduling pre- 
computed sub-models. In this sense, TPWL also fits in a gain scheduling framework, and more specifically a linear 
parameter varying (LPV) approach^^\ which approximates nonlinear system by assembling locally linear models 
obtained by Jacobian linearization at operation points^^l Since LPV method depends on local linear model or first 
order Taylor expansion, therefore, the variation from operation points cannot be very large^^l Other nonlinear 
modeling methods have also been proposed, such as the Discrete Empirical Interpolation Method (DEIM) by 
Chaturantabut and Sorensen^*®^ and the Gauss-Newton with approximated tensors (GNAT) method by Carlberg et 
ak'^^l Both DEIM and GNAT apply POD basis for model reduction. 

Despite successes of POD repotted in the literature, some difficulties and limitations are also present for 
applications to complicated problems in CFD. And we would approach nonlinear modeling with a different idea 
guided by the following principal. We shall retain accuracy of modeling (preferably not losing resolution of the 
original system by reducing the degree of freedom) and achieve efficiency so that the resulting model will be valid 
and economic for a “reasonably wide” variation of parameter’s types and values. 

Previously, we presented a nonlinear reduced model via recurrent artificial neural network training'^*^^ by 
assembling a number of weighted neurons associated with different inputs, in which the weights are based on the 
radial basis functions'^'^l This thinking follows the practice used in metamodeling. The present work utilizes a 
similar concept for constructing a nonlinear model by collecting a set of submodels that are each valid in the 
neighborhood of sampling points. The idea of TPWL comes handy so that we can simply replace the neurons with 
the piecewise linear solutions valid around the sampling points. It is noted that our approach does not belong in the 
realm of model reduction because we do not employ the reduced basis matrix <j> and the dimension in our 
approximate system remains the same as the full order model. 

In the next section, we shall give a detailed description of the mathematical foundation of our method for 
nonlinear modeling that seeks to realize the three properties stated in the beginning. Specifically, we show how the 
submodel at a sampling point is constructed through Taylor series expansion and how these submodels are 
combined through weighting functions to form the final model. Our previous publication'^*"'^ has presented the idea 
and formulation utilizing the first-order Taylor series and the method is called weighted piecewise linear (WPL) 
modeling; the mathematics involved is relatively straightforward. The present paper focus on the second-order 
Taylor expansion and the method is now called weighted piecewise quadratic (WPQ) modeling; the algebra 
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becomes more complicated resulting from the Hessian matrix and seeking its approximation is the subject of this 
study. 


II. Formulation 


In order to determine the evolution of X in Eq. (1) accurately, it is necessary to be able to evaluate f accurately, 
which in general is a complex nonlinear function of X . A straightforward idea is to base on the local information of 
f using Taylor’s series expansion about Xg , as also standard for integrating differential equations. But the validity 
of the expansion is limited to a small region of Xg along the evolution (trajectory) of X . Instead of utilizing just a 
single state, we shall at will approximate f with a series of p locally accurate solutions valid 

around sampling points Xg . and we can assemble them together by weighting functions vr,. 


f(X,U) - X^,(^’«;^0,P»0,,)f,(X,U;Xg ,,Ug ,) 


/=l 


where 


f.(X,U;Xg .,Ug .)=r°^(Xg .,Ug .)+r‘^(X;Xg ,,Ug .)+r"'(X;Xg ,,Ug ,) 






(4) 


(5) 


+ g''(U;Xg .,Ug .)+g' '(U;Xg .,Ug_,) + ... 

with f denoting the function value and f , n>0, the nth derivative with respect to X , and they all are evaluated 
at the sampling point ( Xg ,Ug j ). For example. 


3f 


f (X;Xg,,Ug .)=Z)/(Xg ,,Ug .)-AX, ( X , U ) = ~ , AX = X ~ X g 

a'f 


( 6 ) 


f(^>(x;Xg .,Ug .) = ^ Ax" ■ D"f(Xg ,,Ug ,) ' Ax , ^ ^ I ( X , u ) 


1 


3x' 


where Z)^f(x,u) is the usual Jacobian matrix of f with respect to X and Z)^f(x,u) is the Hessian matrix. 

Evaluation of the Hessian matrix for fluid equations is complicated and expensive, which is the main reason for it 
being scarcely touched in practice. In our earlier work, we also stopped with the first-order expansion. But there is a 
good reason for including the second-order term, such as for enhancing accuracy and allowing a smaller number of 
sampling points while keeping same accuracy, and for automatically including nonlinearity in the formulation. It is 
the primary motivation for the current study in which we seek to develop an effective strategy to approximate the 
Hessian matrix so that the accuracy of second-order expansion is realized while the computational effort for this 
additional accuracy is minimized. 

Similarly , n>0, is the «th derivative denoting the derivative of f with respect to U , also evaluated at the 


sampling point ( Xg . ,Ug . ), 


J’ Uj 
(1) 


df 


8 (»;Xg,,Ug .)=Df(Xg .,Ug .)-AU, ( X , U ) = ~ , AU = U - U g . 


g('>(u;Xg .,Ug .) = ^ Au" ■ D^f(Xg .,Ugp- Au, ^ ^ f ( X , u ) = 


1 


a'f 


(7) 


2 “ ' ' u - ' - ' 

Finally, the weighting function used in this work is based on the radial basis function (RBF)^’^^ of the Euclid 
distance between the state X and the sample (operating) points Xg ^ , given as follows: 
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w. = e 


( 8 ) 


%> 

where is a representative length scale measuring the spread of the solution of the ith sampling point. We set 
<5,= 1 /a/ 50 for all calculations in this study; this choice seems to work well broadly for the variety of problems. No 
attempt has been made to find a better specification of this parameter. It is reminded that the weighting function is 

normalized so that ^ W. = 1. 

In our previous work, we have included the first-order expansions f and g*'* , hence called piecewise linear 
(WPL) in Rewienski and White^*^ ; its power and versatility for the Euler equations for several unsteady 
aerodynamic problems with different input functions are documented'^*"''. Since this model is nonlinear via the 
weighting function, we also call it weighted piecewise linear (WPL) modeling, to avoid the impression for a linear 
model. 

In our current work, we begin by including the second-order expansions, hence the new formulation will be 
appropriately named weighted piecewise quadratic (WPQ) modeling. Including the Hessian matrices Z)^f and Z)^f 

however leads to a formulation with algebraic complication and computational expenses for a fluid dynamics 
systems, despite advantages of increased accuracy and range of trust region. In what follows we shall further explore 
the Hessian terms and attempt to approximate it so that there is no explicit appearance of it, thus significantly 
reducing the complexity of the method and computational cost. 

For illustration, it is sufficient to consider the vector f(x,u) component by component, namely taking j] 




individually. Thus, we shall take the first component for illustration. Also, it suffices to consider the Taylor 
series expansion at an arbitrary sampling point since the formula remains the same for all sampling points. Hence, 
we shall drop the subscript associated with a specific point. And we shall focus on evaluating , while 

the same procedure can be applied to gj^’(x;x^,u^) . Here we have 


1 


y;' '(x;Xo,Uo) = -Ax D /;(x^,Uq)Ax, Ax = x-x„ 


( 9 ) 


where 




f f 

1, -If 1^2 


f f 

^yX2Xj 

•• f 

\,X2X^ 

f f 



f 

dx.dx 

I J 


( 10 ) 


A key observation that renders simplifying the evaluation of the Hessian is explained below. Let us denote 


Ax^z9fy;(Xo,Uo) = [Ji,J2, •••,<] 


( 11 ) 


where a typical element .(x;x.) in the above vector is given by 


J.(x;Xo,u,) = XAt/ (x„,uj 


: Ax. 




* dx, dx 

1 j 


^ dx., dx . 

2 ] 


dx dx 

n J 


Approximate the Hessian derivatives by forward differencing. 


( 12 ) 
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AXj = 

ay; 



dx dx 

dx . 


dx . 

1 7 


y 

X|=x,+Ax,,x^.(yjil) 

J 

dV, 

At — 

a/: 


ay; 

3x, dx . 

ZXV 2 

dx 


dx. 



J 

X 2 =^2 + Ax2 a - ( y 9^2 ) 

J 


Ax = 

a/: 


J_A 

dx dx 

n 



dx . 

n y 


y 

~^n ) 

J 


By substitution of the above approximations, d . becomes 


J (x;x,,u„) = (x„,u„) 


+ Ax,(^2+^2’W2)’«o)-4„(Xo’«o) 


+ ••• + 


= Z /l,x . ’ Uo ) - «/l,x , (*0 ’ «0 ) 

/=1 

Define a mid-point, X , 


- 1 A x 

X = — (x + x„) = x„ H x = x„ + Ax„ 

2 ° 2 


0 0 


( 13 ) 


(14) 

(15) 


The first term in Eq. (13) is further expanded at this midpoint, it is then contracted to give a considerable 
simplification. 


X/ 

d . = n—^ 

^ dx. 

j 


dx 


(16) 


We call attention to the following: (1) only the first derivatives, no longer second derivatives, are needed, and 
(2) they are evaluated only at two states, (x,Xg) . 

Finally, we have the second-order expansion now in terms of a difference of first derivatives respectively 
evaluated at the mid-point and base states 


= ^Ax^dVj(x„,u„)Ax = 


fl 

~2^ 


-(x,u„)--(x,,u„) 


(17) 


Ax. 


Replacing the subscript “1” with “j”, we get the expression for the yth component 


y;^'^(x;xo,Uo) = |X 

^ /=1 


5 / _ 5 / 


Ax. 


(18) 


Thus, the expression for the full vector f ^^'(x;x. ,U„) is now obtained. Remarkably, these first derivatives are 


already available when evaluating the first-order Taylor terms or f^‘^(x;x„,u„). This amounts to using the 
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same information required for the first-order Taylor expansion, except the Jacobian matrices are now evaluated at 
the average state, instead of the sample points. Similarly, can be obtained easily. Thus, the yth 

component is summarized as follows. 

//x;x„,u„) = //x„,Uo) + (l-^)X^(Xo,Uo)Ax,. + ^X^(X’«o)^,' 


l — i , l — l , 


(19) 


Now it is instructive to escape the above abstraction and lengthy algebraic manipulation of the n-dimensional 
system and consider a simpler one for bringing some clarity. Let us consider the following two-equation system. 


a 

X 



dt 

y 


f2i^,y,u) 


As before, the local Taylor’s series expansions about a sampling point (Xg ; ,To , , ) reads 

/i 


f(x,y,u)- 


A 


(20) 


(T,T,M)-f^°’(Tg.,3;g,,.,Mg,,.) + f^‘’(T,3;;Tg.,3;g,,.,Mg,,.) + f^"’(T,3;;Tg.,3;g,,.,Mg,,.) 

+ g^‘’(M;To.,To,pMo,) + g^"’(M;Tg.,3;o,.,Mo,) + --- (21) 


where 


and 


f^'’(X,T;Xg.,Tg,pMg,,.) = 


1 

1 


Ax, 


Ax, 



1 

>> 

1 



5 

. . 


1 

1 

O 

1 


f^"’(T,T;Tg_,.,To,,-,Mo,,) = ^ 


Ax, At, 


Ax, At, 


/i,x. /i,: 

f f 

J \^xy ^-.yy 


f 2 ,xx f 2 ,yx 

•f2,xy ^2,yy 


“1 

Ax, 


4(Xo,,Jo.i."o,') 

. . 



Ax, 



. . 



With t - 1 , 2 , etc. 

ox oxox 

Similarly, 


g^‘’(M;Xo_,,To,,-,Mg,,)^ 


fl,u 

fl,u 


Am,, Am, = M-Mg, 


(X„,, Jo 


g^"’(M;Xo,,To,-,Mo,) = ^ 


/i,„, 

fl,u 


(Am,)" 




The second-order Taylor term of can be expressed as (from Eq. (23)) 


( 22 ) 


(23) 


(24) 


(25) 
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( 26 ) 




where 


J, =Ax,./i 
fi?2 = Ax,. /i 


2 

I +Ay.f. I 

^ +^yifuyy\ 


(27) 


Applying forward differencing now to approximate the second derivatives, we get 

= )+fu(y’\p \i )] 

(28) 

d2 = [A,y(x-,yo,,u,,)-2f^^^ (Xq . , Jo . ,u^.) + f^/y; )] 

, we get the following trapezoidal approximation for and , as 




Expanding the derivatives about an average state 
follows. 


(29) 


and 


(30) 


A = [2/,, (x , J , Mo )-2f^^ (Xo , Jo,, , Wo, )] 

A = [2/;,,(A, J,Mo,)- 24 ^(Xo,, Jo,,Mo,)] 

The formulas for ,gj and g 2 can be obtained similarly; it is even simpler if only one variable u is considered 
in this paper. Putting both the first and second Taylor terms together, we get an approximation of at the ith 
sampling point (^o ,, Jo,,Mo,) 

/l(x,J,M) = /i(Xo,,Jo,,-,Wo,) + /i,,(x,J,Mo,)Ax,. + /io,(x,J,Mo,)Aj,. + 4„(Xo,,Jo,i,w)AM,.+... (31) 

^re 


where 


^ = y = \(y+yo,A u = ^{u+u^.) 

is obtained similarly, by simply replacing the subscript “1” in with “2” 


(32) 


The formula for -r - -r 

Then, the complete WPQ for the system Eq. (20) reads 


dx 

¥ 


p 

f(x,u)-^tv.(x,u;Xo,,Uo,) 

/=1 


A 

/2 


(x,u;x u ) 


(33) 


where X = (x, j)^ and U is a scalar function if the formula in Eq. (31) is to be used directly. This is all what 
the WPQ modeling is about and it is ready for a usual time integration provided an input data is given. This is the 
basis of modeling where a baseline (arbitrary) input function is employed and the integrated is performed; the 
process is called training. The important piece that shows its power — main impetus of our work, lies in the 
application of this modeling (training) to a different kind of input functions without having to recalculate the 
original f for each new input, which is an expensive proposition. 


III. Computational Results 

We will show in this section the results from using the modeling method detailed above, for problems involving 
a scalar equation, a system of two equations, and Euler equations with an input that differs vastly in form and 
magnitude from the baseline. The purpose is to show the power and applicability of the method, with regards to 
its accuracy, cost savings, and versatility/robustness. As a benchmark solution for comparison for its validity, 
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integration of the original system will be treated as the “exact” (albeit numerical) solution, which is obtained by 
using a four-stage Runge-Kutta method. 


A. Scalar Equation 

We employ a nonlinear scalar equation with an input function u to facilitate understanding of the WPQ modeling. 

— = f{x{t),u{t)) = 5smx + x + u (34) 

at 

The WPQ modeling solves 


dx 

ot (35) 


a. = 5cosx. + l, b. = \, A =-2.5sinx 

l l ^ l ^ l I 


Here the exact Hessian is used in Eq. (35) as it is readily available; its results also serve as bases for comparison 
with that of approximate Hessian described above. Specifically, the approximate Hessian developed in previous 
section gives 


dx 

■^ = Z [/(^, ’ -Xi) + b^(u- u . )] 
a. = 5cosx. + 1, x = (x + x)/2 

For illustration, we arbitrarily choose the following input for training 

m( 0 = sin6^ + cos6^ 


(36) 


(37) 


Then the original equation, Eq. (34), is integrated over a specified interval, say (0,T=27t), providing a solution at 
N discrete point (state) and the coefficients, (a.,b.,h.) and (a.,b.). Out of these N points, we pick p sampling 

points with their coefficients. The sampling may be arbitrary, we chose an even distribution of points in (0,7); other 
sampling strategy and its effect are presented in our previous study'^*'*l These trained coefficients now can be used 
in Eq. (35) or (36) with the specified sampling points and the integration can proceed as usual to give the 
approximate solution. 

First, we shall use the exact Hessian, Eq. (35). Figure 1 shows the solutions using various sampling points and 
the comparison with the “exact” solution. The solution with only one sampling point explodes quickly, but with 
additional sampling points it converges to the exact solution well. 



Figure 1. Comparison of various WPQ solutions with the “exact” solution. 

Next we use the approximate Hessian and the solution is shown in Figure 2, here the solution with even just one 
sampling point surprisingly follows the correct trend (compare with Figure 1), the approximate Hessian gives a 
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faster convergence to the exact solution. In Table 1, one can clearly see a drastic improvement in accuracy by 
including the second-order expansion term, and the approximate Hessian yields more accurate solutions than the 
exact Hessian. Moreover, all three methods shows convergency as more submodels (sampling points) are used, i.e., 
as p is increased. The error is measured, throughout this work, by the following formula. 


Error 



\ : exact solution 


(38) 



Figure 2. Solutions by the approximate Hessian using different numbers of WPQ submodels. 


Table 1. Comparison between first-order and seeond-order Taylor expansions. 


PW Elements 

First-order 

Exact Hessian 

Approximate Hessian 

1 

NaN 

NaN 

8.37e-2 

2 

2.27e-2 

1.92e-2 

6.15e-3 

5 

1.20e-2 

1.78e-2 

5.50e-3 

10 

8.05e-3 

6.82e-3 

1.67e-3 

20 

4.19e-3 

8.21e-4 

1.74e-4 

50 

1.15e-3 

5.15e-5 

1.13e-5 

100 

3.91e-4 

7.16e-6 

1.61e-6 

150 

2.52e-4 

2.95e-6 

6.87e-7 

300 

1.70e-4 

1.28e-6 

3.12e-7 


Next, we use the training data obtained with Eq. (37) and extend it to a new input function, 

M(t) = 5sin^6t (39) 

The result with 10 WPQ submodels gives an excellent comparison with the exact solution, as shown in Figure 3; the 
training solution is also included for reference, showing a significant difference between these two solutions, thus 
confirming the versatility of the approach. 
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Figure 3. Results of using 10 WPQ submodels for input Eq. (38). 

Next, we consider a more complicated input function, which can produce an exploding increase in solution x. 

u{t) = sin^ 6t (40) 

If jJ. is small, a sudden increase in the solution can occur, as seen in Figure 4 for p.=0.25. The solutions by the exact 
and approximate Hessians are again included for comparison, in which 1 WPQ again already shows main dynamics 
of the solution and the 5 WPQ solution essentially coincides with that of 10 WPQ. Noticing again the significant 
effect of the new input function, compared to the baseline solution. This example shows the reliability of the method 
for a case with an increasing amplitude in the solution. 



Figure 4. Solutions by tbe approximate Hessian using different numbers of WPQ submodels, tbe training and 
exact solutions are included for comparison. 


B. System of Two Equations 

Here we specify the 2-equation system considered previously 


a 

X 


fM,y,u) 

dt 

y 




with 


( 41 ) 
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( 42 ) 


/j(x,j,M) = sin(x + j) + M 

f^{x,y,u) = sm{xy)-u 

The input function for training is same as Eq. (37). 

Using the trained data, we now apply the model to a new input 

u{t) = sm^6t (43) 

The solutions by the first-order and second-order Taylor expansions with 10 submodels are compared with the exact 
solution for the input Eq. (43) in Figure 5; the first-order solution has a significant discrepancy (even though giving 
the correct trend) in the x component from the exact solution, but the second-order expansion provides excellent 
accuracy for both x and y. We note that the first-order solution does converge to the correct solution, by using more 
submodels, as shown in the same figure. 




Figure 5. Left: eomparison of solutions by the first-order Taylor (WPL) and second-order Taylor (WPQ) 
expansions with the exact solution, using 10 submodels. Right: convergence of WPL using 5, 10 and 40 
submodels. 

Using yet another input given in Eq. (40) for the same training data, we get a vastly different result, as shown 
in Figure 6, with training data included for reference (please note that the vertical coordinate is adjusted in both 
plots for clarity). In the input with exponential function, the solution again shows increasing amplitudes in both x 
and y components; the 10 WPQ model gives an excellent comparison with the exact solution, even for cases far 
from the baseline. 




Figure 6. Solutions using 10 WPQ submodels with new inputs, with the training solution (black and red 
curves) included for comparison. Left: input Eq. (43) and Right: input Eq. (40). 

C. Euler Equations 

Let us consider the 2D Euler Equations written as: 
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(44) 


= R(q,u,v) 

where V is the volume of a computational cell, q consists of the conservative flow variables, R(q,u,v) is the 
residual representing spatial derivatives of flow flux functions for a set of input conditions, u and v, which may be 
the structural displacement and velocity, such as in a pitching/plunging motion of an airfoil or the flow condition 
such as the Mach number. The spatial discretization contained in the residual term R(q,u,v) can be determined by a 
state-of-the-art numerical flux method; here we employ the AUSM^-up scheme'^*^^ to approximate the flux function 
and the van Albada limiter in the MUSCL interpolation of primitive variables. 

At the time of finishing up this paper, we have not been able to get the WPQ implemented in the code to work 
properly. Hence we shall show in the following some typical results presented in the past publication'^*"'', which were 
obtained by using the weighted piecewise linear (WPL) model, in order to illustrate the efficacy of the approach for 
a fluid dynamics system. The interested reader is encouraged to refer to our previous paper for details. The WPL 
model is expressed by 



djvq) 

dt 


= X + A,(q - q,) + B; ,(u - u,) + - v .)} 


(45) 


where the Jacobian matrices are 


A=^ 

' 3q, 


B. = - 

3u 


B 

3v 


Q = [q,u,vf 


(46) 


The right-hand term of Eq. (45) gives a weighted combination of p locally linearized models. The evolution of this 
system is now completely a function of the baseline (training) solutions at the sampling points; it is computationally 
tractable, requiring only the matrix-vector multiplications and the evaluation of the weights. The evolution of q 


now completely bypasses flux evaluations of the entire flow domain as needed in the full Euler equations for each 
time step. 

We shall see some unsteady aerodynamic characteristics associated with some standard test cases in the 
AGARD report'*"*', specifically the so-called “CT2” and “CT5” cases. These two cases have different flight Mach 
numbers and pitching parameters, as given in Table 2. The pitching motion of the airfoil is described by its angle of 


attack varying in time from the mean value : 


a{t) = a +a.sm(ot = a +a„smlkt, 

^ m 0 m 0 ’ 




(47) 


where (X^ is the oscillation amplitude, k is the reduced frequency defined by the free-stream velocity [/« and the 
chord length c. The pitching motion is centered at Xm- 


Table 2. Parameters of CT2 and CT5 eases 


Case 


(Xm 

ao 

k 

Xm 

CT2 

0.6 

3.16° 

4.59 

0.0811 

0.273 

CT5 

0.755 

0.016° 

2.51 

0.0814 

0.25 


First, we validated the accuracy of the CFD procedure against the experimental data'*"*'; the comparison of 
aerodynamic forces, in terms of lift and moment coefficients, in a cycle is given in Figure 7 for both the CT5 and 
CT2 cases with a sequence of grid refinements, showing the CFD results in good agreement with the data. More 
importantly, both test cases exhibit a strong nonlinear behavior in the moment coefficient and nonsymmetrical 
forces between downward and upward motions. It is noted that CT2 has a significantly larger excursion of about 10° 
in pitching than CT5 with 5°, thus it is expected to show stronger nonlinear effects. In contrast, CT5 operates at a 
higher Mach number. Both have a varying supersonic pocket in the flowfield during pitching, also manifesting 
nonlinear effects. 
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a 


(a) CT5 


o 



O 


(b) CT2 

Figure 7. Validation of CFD solutions for the AGARD CT5 and CT2 problems. 


In what follows, we separately discuss the results from the nonlinear modeling for the CT5 and CT2 problems, as 
they exhibit rather different flow characteristics even subject to the same input parameters. 


Cl. CT5 

After the training at the baseline condition ao= 2.51°, k = 0.0814, the so-built WPL model can be used for other 
conditions. Shown in Figure 8 are the results of the full CFD model and the nonlinear model using 31, 21, and 1 1 
WPL solutions with «o= 3.5°, k= 0.0814; it shows that both 21 and 11 WPL models can also faithfully capture 
nonlinear features of the original CFD system. In contrast, the linear results, obtained by retaining only the zeroth- 
order Taylor expansion around the steady operation point, completely miss the nonlinear phenomenon in the original 
system. As the linear model does not contain any time varying information in the model other than via input 
function, it is clear that collecting time-varying submodels into a complete model is essential and is a key to 
maintaining accuracy of the model. After the training, the so-built WPL-based model can be used for other 
conditions. 
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■ Euler 


O 



Figure 8. WPL results eomputed for a different pitehing motion «o = 3.5°, k = 0.0814 from the baseline 
(training) data of «o = 2.51° and k = 0.0814. 

Now we consider the scenario of having a different flight motion from the baseline. In this case we impose a 
new input signal but with a cubic sinusoidal function 

a-a^+a^sm^{2kt) (48) 

with the reduced frequency k and ttm unchanged, but with different amplitudes ao=3.5°, 1.0°and 5.0°. It is noted 
that the cubic sinusoidal input is expected to give rise to a more complicated flight motion. Also, the amplitude of 
5.0° is significantly larger than the baseline value of 2.51°. 

Figure 9 shows that the WPL-based model still hold its accuracy even the flight motion is both qualitatively 
(cubic function) and quantitatively (amplitude) different from the training motion. 

Euler 





■ Euler 



(b) Uo = 5.0° 
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Figure 9, WPL results with Eq. (48) for amplitude Uo = 1.0° (a) and and 5.0° (b) departing from the baseline 
ao = 2.51° 


The proposed modeling, as mentioned earlier, maintains the same degrees of freedom as the full model, unlike 
the model reduction methods, it is interesting to compare the detailed flowfield predicted by the model with the 
solution calculated by integrating the full Euler equations. Figure 10 shows pressure contours of both solutions side 
by side; they are nearly indistinguishable. 



Figure 10. Flowfield expressed by pressure contours at a time slice (t=T/4) with a cosine signal (no = 3.5°, At = 
T/60) (Left: Euler, Right: WPL). 


Finally, we compare in Table 3 the computational cost for performing a full-order CFD analysis and the 
additional work required to construct the WPL-based nonlinear model. One should note that this extra overhead is 
one time only; any computation for other conditions will require a minimal additional effort to perform matrix- 
vector multiplications and time integrations — only 1/50 (for the current setup) of the cost incurred by a full CFD 
solution is needed for another new condition by the model. Hence, in this sense, the so-constructed nonlinear model 
is not only accurate but also computationally economical for parametric studies and design optimization. The more 
calculations performed, the more savings realized. 


Table 3. Comparison of CPU times by full model, WPL model construetion, and running a eonstrueted WPL 


model on an Intel 17 CPU. 

Model Choice 

CPU Time (min) 

Full Model (Euler) 

15.03/period 

WPL Model Construction 

8.533 (for 60 WPL samples/period) 

WPL Model 

0.324/period 


C2. CT2 

For another pitching airfoil problem, CT2 we further apply the previously established WPL-based nonlinear 
model. This problem is more complex aerodynamically than CT5 in the sense that the aerodynamic forces display 
stronger nonlinearity, as seen in Figure 11, where the lift force has a noticeable nonlinear behavior near high angles 
of attack, exhibiting asymmetry with respect to angles of attack. 

In Figure 11, we show the performance of the WPL-based solutions against the full Euler solutions of different 
pitching amplitudes. At a small amplitude («o = 1 °), the flow behaves linearly and the linear model is close to the 
Euler solution. However, when the amplitude is moderate, nonlinearity first appears in the moment coefficient; then 
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the lift coefficient exhibits a strong nonlinear phenomenon with a double loop, suggesting that a second, higher- 
frequency flow structure is embedded in the main flow. The linear model clearly is incapable of predicting the 
correct flow behavior. 




Figure 11. Comparison of Euler, nonlinear WPL-based model and linear solutions of the CT2 airfoil at 
different pitching amplitudes, with no = 8°. 

A snapshot of the detailed flowfield is displayed in Figure 12 for the training condition («o = 8°); no discemable 
differences between the Euler and WPL-based solutions can be found. 




X 


P 

2.40 
2.20 
2.00 
1.80 
1.60 

1.40 
1.20 
1.00 
0.80 
0.60 


Figure 12. Pressure contours of AGARD CT2 at a time instant, r/4, for the airfoil pitching at an = 8° (Left: 
Euler, Right: WPL). 


Again, we are interested in seeing how the model so constructed performs outside the training trajectory for the 
CT2 problem. With the previously trained model for the trajectory prescribed by Eq. (47) with ao = 8°, and = 
3.16°, we now apply it to a new motion defined by Eq. (48), which contains high (triple)-frequency contents 
resulting from the cubic function and inflection points at t=0, T/2 and T. 
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Training Input Eq.(45) 



Figure 13. New input trajeetory defined in Eq. (46). 


A vastly different behavior is revealed in both the lift and moment coefficients in Figure 14, showing sharp 
changes at the mean angle «,„= 3.16° in both upward and downward phases. At the knot point (a == 6°) in Q, where 
Cl trajectory crosses, Cm however does not find a particularly unusual behavior except having double values. 




Figure 14. Lift and moment eoefficients responding to Eq. (48) input with Uo= 8°. 


Finally we show the flowfield at time slices t = T and 3772, both at knot points but with opposite pitching 
directions, as shown in Figure 15; the results computed by the model are nearly identical to the CFD results in every 
detail of the flowfield, also indicating a nonsymmetry here even though the airfoil is at the same angle (a = 3.16°) at 
both times. 


17 

American Institute of Aeronautics and Astronautics 



- 0.5 0 0 5 1 - 0.5 0 0,5 1 

X X 


(b)r = 3772 


PFigure 15. Detailed flowfield eomparison by tbe WPL-based model and CFD (Left: Euler, Right: WPL), «o = 
10 °. 


IV. Concluding Remarks 

In many scientific and engineering endeavors, there occurs needs for having a mathematical model that can 
faithfully describe the physical phenomena under conditions of changing parameters. Examples include design 
optimization where some geometrical parameters or flow conditions will be required to be changed in order to arrive 
at an optimal result; similar scenarios also appear when exploring fluid dynamics phenomena resulting from 
different flow parameters or geometrical effects. A reliable and useful model must be sufficiently accurate to capture 
important characteristics of the problem under consideration — generally nonlinear in nature and must also remain 
valid for a sufficiently wide range of parametric changes. 

In this paper, we presented an innovative method for constructing a versatile predictive model and we 
demonstrated for problems for simple nonlinear scalar equation to Euler equations. The method is a nonlinear 
combination of Taylor series expanded about a set of sampling points (states). The model is trained with a baseline 
input function and is capable of predicting solutions for different inputs. Built upon our earlier publication where the 
first-order Taylor expansion was used, we studied the feasibility of including the second-order Taylor terms, which 
involve complex and computationally expensive Hessian matrices for a system of equations. An approximation of 
the Hessian matrix was proposed. The accuracy and robustness have been demonstrated for the scalar equation and a 
two-equation system using this approximate Hessian, and for the Euler equations using the first-order Taylor 
expansion. Work to implement the approximate Hessin in the Euler system is continuing. 
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